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The wetting properties of solid surfaces are significant in oil/gas and liquid displacement processes. It is 
difficult for hydrophobic fluids to permeate channels filled with hydrophilic particles and an aqueous phase, 
and this is thought to be the primary cause of low yields in low permeability reservoir operations. Using 
three-dimensional lattice Boltzmann simulations, we show that particles with hydrophobic and hydrophilic 
patterned surfaces can greatly improve hydrophobic fluid permeation. Specifically, a hydrophobic fluid can 
easily access micro-channels in the hydrophobic regions, which extend rapidly even to the hydrophilic 
regions and accelerate hydrophobic fluid escape. This work enriches understanding of multiphase flow in 
porous media at the pore scale and fracture conductivity and is expected to have great significance in the 
exploitation of low permeability reservoirs and shale gas. 

Low permeability reservoirs are characterized by very low matrix permeability, minute pores and natural non- 
connected fractures, and thus have low natural industrial capacity. To improve the recovery efficiency of oil 
and gas, hydraulic fracturing is usually applied after wells are drilled 1 . Hydraulic fracturing is implemented by 
pumping large volumes of an aqueous phase containing suspended proppants into reservoirs under high pres- 
sure. The rocks around wells are split or fractured and the proppant particles infiltrate and support the fractures. 
Usually, the surfaces of proppant particles are hydrophilic so they can penetrate the fractures under the aqueous 
phase flow. Because of the strong attractive interactions between the aqueous phase and particles, the proppants 
will trap the injected aqueous phase and then form barriers preventing outflow (diffusion) of hydrophobic oil/gas. 
This phenomenon is called water blocking damage 2 and is the primary cause of low yields in low permeability 
reservoir developments during the middle-late stage 3 " 5 . It is clear that a key objective is to keep the proppant - 
water interaction strong enough, while allowing the hydrophobic fluid to flow (diffuse) rapidly near the proppant. 
Through the development of new synthesis techniques 6 ' 7 , it is becoming easier to decorate particle surfaces 8 . We 
propose that the presence of hydrophobic patterns on hydrophilic particle surfaces may greatly enhance per- 
meation of the hydrophobic oil/gas in micro-channels. Because there would still be sufficient hydrophilic regions 
on the particles, interactions between the particles and aqueous phase would still be sufficient to enable its 
penetration into a fracture. On the other hand, there may be inhomogeneous wetting properties on the surfaces 
of natural rocks. Thus, it is also of great importance to study oil/gas permeation through fractures in such rocks. 

In this paper, using a three-dimensional (3D) model, in which the channel is packed by particles with 
hydrophobic and hydrophilic patterned surfaces, we show that wettability-patterned particles can greatly 
enhance hydrophobic fluid permeation through a channel. We used a lattice Boltzmann method to perform 
numerical simulations and found that a hydrophobic fluid could easily displace an aqueous phase from the 
hydrophobic regions and opened micro -channels. In addition to hydrophobic fluid infiltration, these micro- 
channels extended rapidly even into the hydrophilic regions, resulting in elimination of blocking by the aqueous 
phase in the particle-packing region. 

The geometry of the system is shown in Figure 1(a), together with a snapshot of some aqueous phase fluid. The 
length of the channel was 480 lattice units, while its cross-section was a 120 X 120 square array. Arrayed in two 
layers, 18 spherical particles were closely packed from x = 130 to 210 in the channel. The diameter of each 
spherical particle was 40 lattice units, and thus the volume of the packed porous medium was 80 X 120 X 120. The 
non-slip bounce-back scheme was adopted for the particle boundary and channel wall, and a pressure boundary 
condition was employed at the inlet and outlet to drive fluid flow. Initially, the aqueous phase filled the space from 
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Figure 1 | (a) Schematic diagram illustrating water blocking, in which the light yellow denotes the aqueous phase in the initial equilibrium state, (b) 
Sketch of the wettability-patterned particles with two hydrophilic and two hydrophobic regions packed in a channel, (c) Sketch of wettability-patterned 
particles with four hydrophilic and four hydrophobic regions packed in a channel, (d) Sketch of the mixed hydrophobic and hydrophilic particles 
packed in the channel. The red parts indicate hydrophobic regions, with an intrinsic contact angle of 117.7 degrees, whereas the blue parts indicate 
hydrophilic parts, with an intrinsic contact angle of 32.8 degrees. 



100 to 240 in the channel, while the hydrophobic fluid filled the 
remaining space, and the same pressure was assigned at the inlet 
and outlet. The system was initially equilibrated during (typically) 
10,000 time steps. We computed the velocity of fluids through a 
cross-section (x = 60), as shown in Figure SI in the supplementary 
information. We can see that after 8000 time steps, the velocity 
evolved to achieve equilibrium. We therefore conclude that the sys- 
tem achieved an equilibrium state after 8000 time steps. We show a 
snapshot of fluid distributions in the channel in Figure 1 (a) at the 
10,000th time step, after which the process of fluid flowback was 
simulated and evolution of the fluid flows in the channel was 
recorded as a function of time by decreasing the pressure at the outlet. 

We studied the influence of particle surface wettability on the 
fluids flowing in the channel. Two well-defined types of wettability- 
patterned surface are shown in Figure 1 (b) and (c) and we performed 
simulations with mixtures of hydrophilic and hydrophobic particles. 
The mixed forms are shown in Figure 1 (d). For comparison, we also 
performed simulations with fully hydrophobic and fully hydrophilic 
particles. In the simulations, the intrinsic contact angle of the hydro- 
philic surface was about 32.8 degrees and the intrinsic contact angle 
of the hydrophobic surface was about 1 17.7 degrees. The channel wall 
was hydrophilic and had a contact angle of about 11.4 degrees. 
(Simulation details are shown in the supporting information.) 

Results 

When we reduced the pressure at the outlet, fluids flowback began. 
The evolution of the hydrophobic fluid flow through the cross-sec- 
tion (x = 60) for five cases was recorded and the velocity as a function 
of time is shown in Figure 2. It can be seen that the velocities first 
increase and then decrease to various extents. Fluctuations can 
clearly be seen in the patterned- wettability cases and the mixed case, 
where the hydrophobic fluid drives the aqueous phase away from the 
porous medium. The velocities each finally achieve steady state after 
100,000 time steps. However, for hydrophilic particles, the hydro- 
phobic fluid velocity reduces to a very small value, meaning that the 
hydrophobic fluid is blocked by the aqueous phase. 

Because particles with wettability-patterned surfaces significantly 
reduce the interaction between the aqueous phase and the particle 
surface, the aqueous phase in the hydrophobic region can be readily 
cleaned up. Our simulations confirm that a patterned surface can 
markedly enhance the ability of a hydrophobic fluid to break through 
water blocking in micro-scale pores. 

Discussion 

We have calculated the aqueous phase volume fraction in the pore 
space of the particle -packing region (denoted by S w ), namely from x 



= 130 to 210 lattice units in the channel. Figure 3 shows S w changing 
with evolution time. 

In all five cases, S w began to decline after we reduced pressure at 
the outlet, and tended to be stable after about 50,000 time steps. The 
stable volume fractions for the five cases clearly differ. S w for the 
hydrophilic particles was remarkably higher than for the other four 
cases, meaning that there was more aqueous phase entrapped within 
the porous region. Because the hydrophilic regions of wettability- 
patterned particles retained a small quantity of aqueous phase adja- 
cent to the channel wall, S w was in the medium range. Conversely, the 
fully hydrophobic particle eliminated most of the aqueous phase and 
displayed the greatest conductivity. This explains why the fully 
hydrophobic particle shows the largest hydrophobic fluid velocity 
in Figure 2. Meanwhile, the volume fractions of the wettability-pat- 
terned particles in Figure 1(b) and 1(c) and the mixed case in 
Figure 1(d) are almost the same, so that the three cases had similar 
conductivity and display the same fluid velocity in Figure 2. Images 
illustrating the analyses above can be found in Figure 4. 

To display the distribution of entrapped aqueous phase in the 
particle-packing medium, snapshots of a cross-section along the x 
axis at)/ = 60 at different times (0, 10,000, 50,000, 100,000, 150,000, 
170,000) are shown for all five cases in Figure 4. All the systems were 
in the equilibrium state at t = 0. It is clear that some aqueous phase 
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Figure 2 | Comparison of the velocity of hydrophobic fluid through a 
cross-section (x = 60) for different cases. Black solid squares ( ■ ), red 
solid circles (•), blue triangles (A), green inverted triangles (T) and 
magenta left triangles (M) correspond to hydrophobic particles, particles 
modified as in Figure 1 (b), particles modified as in Figure 1 (c), mixed 
hydrophobic and hydrophilic particles, and hydrophilic particles, 
respectively. 



SCIENTIFIC REPORTS | 4:5738 | DOI: 1 0.1 038/srep05738 



2 



I 1 I 1 I 1 I 




0 50000 100000 150000 



Time step 

Figure 3 | Volume fraction in the particle-packing region changing with 
evolution time. Black solid squares ( ■), red solid circles (•), blue 
triangles (A), green inverted triangles ( ) and magenta left triangles (A) 
correspond to hydrophobic particles, particles modified as in Figure 1 (b), 
particles modified as in Figure 1 (c), mixed hydrophobic and hydrophilic 
particles, and hydrophilic particles, respectively. 

near the particles formed water blocking in each case. With system 
evolution, the hydrophobic fluids flowed towards the particle-pack- 
ing regions. At the 10,000th time step, we can see that the hydro- 
phobic fluid had driven the aqueous phase inside the micro-scale 
pores in all cases. At the 50,000th time step, the hydrophobic fluid 
had broken through the particle -packing medium for particles with 
hydrophobic surfaces (a), patterned surfaces (b) and mixed hydro- 
phobic and hydrophilic particles (d). At the 100,000th time step, for 
particles with hydrophobic surfaces (a), patterned surfaces (b) and 
(c), and mixed hydrophobic and hydrophilic particles (d), the hydro- 
phobic fluids drove the aqueous phase out. For case (e), the hydro- 
phobic fluid was blocked by the aqueous phase near the particles. At 
the 150,000th time step, in both cases (a) and (b), the hydrophobic 
fluids broke through the water blocking and flowed out of the chan- 
nel. At the 170,000th time step, in both cases (c) and (d), the hydro- 
phobic fluids also broke through the water blocking. However, in 
case (e), where hydrophilic particles were used, the hydrophobic fluid 
could not break through. In cases (b), (c) and (d), the hydrophobic 
regions were the same as the hydrophilic regions on the particle 
surfaces, whereas the times that hydrophobic fluid broke through 
water blocking differed. This indicates that the particles in 
Figure 1(b) were more effective with respect to hydrophobic fluid 
breakthrough. To understand the distribution of immiscible fluids in 
the pore space of case (b), we have drawn snapshots of the cross- 
section (x = 150) at different times in Figure 5. 



Figure 5 gives the details of distribution of immiscible fluids in the 
micro- scale pores for case (b). We can see that the hydrophobic fluid 
started to open up micro -channels in the hydrophobic region of the 
center particles at the 4,200th time step. The micro -channels grew 
rapidly and covered almost all hydrophobic regions of the center 
particle along with the hydrophobic fluid penetration. At the 
8,000th time step, new micro- channels began to form on the hydro- 
phobic region of the outermost particles. Meanwhile, the middle 
micro- channels extended outwards and after uniting to form vertical 
micro- channels at the 8,200th time step, they gradually broadened, 
even extending into the hydrophilic regions, and spread into the 
spaces near the channel wall at the 10,000th time step. The micro- 
channels achieved steady state at about the 100,000th time step. In 
the hydrophobic region, the interaction between the aqueous phase 
and solid wall becomes weak and the hydrophobic fluid can easily 
remove the aqueous phase from hydrophobic regions, until the 
hydrophobic fluid breaks through the water blocking and flow out 
of the channel (shown in Figure 4). 

In summary, using a 3D model including a channel and particles 
having surfaces with different wetting behaviors, we show that, in 
contrast to hydrophilic particles, wettability-patterned particles 
weaken the attractive forces between the aqueous phase and the 
particle surface and enhance channel conductivity for the hydro- 
phobic fluid. For hydraulic fracturing used in the recovery of oil 
and gas, use of proppant particles with hydrophobic and hydrophilic 
patterned surfaces, is expected to eliminate water blocking. This 
work enriches our understanding of multiphase flows in porous 
media at the pore scale and fracture conductivity. It is very important 
to study oil/gas permeation in rock fractures, and is significant for the 
exploitation of low permeability reservoirs and shale gas. 

Methods 

The lattice Boltzmann method is well- suited to simulate the movement of fluids 
through porous networks 9 " 18 . We used the Shan-Chen multi-component multi-phase 
model 19,20 to study gas/oil flow through an idealized porous medium in micro-scale 
pores. This model involves a single relaxation time in the Bhatnagar-Gross-Krook 
(BGK) collision operator 21 . The time evolution of the o th component in this model 
can be written as 

f?(x+c i At,t + At)-f?(x,t)= - 1 ]ff{x,t)-fi"{xjt)\, (1) 

where At is the time step,/? (x,t) is the density distribution function of the o th 
component in direction c ; at x, f° ,eq (x,t) is the equilibrium distribution function of 
the o th component. x a is the relaxation time of the o th component characterizing the 
collision processes by which the distribution functions relax towards their equilib- 
rium distributions. 

In the three-dimensional nineteen-velocity model, the equilibrium distribution 
function,/' 7 '^*,^), depends only on local density and velocity and can be chosen as 
having following form: 



(a) (b) (c) (d) (e) 




Figure 4 | Snapshots of the fluid distributions for differing particle wettability and evolution time, (a) The case with fully hydrophobic particles, (b) The 
case with wettability-patterned particles as in Figure lb. (c) The case with wettability-patterned particles as in Figure lc. (d) The case with 
hydrophobic and hydrophilic particles mixed as in Figure Id. (e) The case with fully hydrophilic particles. 
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Figure 5 | Snapshots of a cross-section (x = 150) at different times for particles with two hydrophilic and two hydrophobic regions as in Figure lb. 
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where the discrete velocities are 

!(0,0,0)c 1 = 0 
( + l,0,0)c,(0, + l,0)c,(0,0,±l)c i=l,. ..,6 
(±l, + l,0)c,( + l,0,±l)c,(0, + l,±l)c z = 7,...,18 

and the weight factors are 

r 1/3 ;=o 

C0i=< 1/18 i=l,2, ••-,6 . 

[ 1/36 f = 7,8, • • • ,18 



(2) 



(3) 



(4) 



The mass density p a and momentum density p a u a of each component fluid are 
calculated for each node via the following equations 



The relaxation time tunes the kinematic viscosity of the a th component by 

v" = (r"--)c 2 s At. 



(6) 



(7) 



The interactions between gas, liquid and rock control the flow and distribution of 
fluids in porous media 22 " 26 . In the Shan-Chen multiphase model, the fluid-fluid 
interaction, which is proportional to the product of a function of the particle number 
densities, can be written as 



Ff(x) = - G c p a (x,t) Wip a (x + CiAt,t)ci, 



(8) 



where F f(x) is the fluid-fluid interaction force, and the strength of two-phase 
interaction is denoted by G c , which is negative and represents the attraction between 
the two phases, o and o denote two different fluid components, respectively. w f is a 
weight factor, depending on the velocity model. It is simple to describe the interaction 
between a fluid and a wall by introducing an extra interaction force, and this method 
is first used by Martys and Chen 27 . The idea is to create an analogue to the particle- 
particle interaction force used to induce phase separation, and it is described as 

Kds^) = ~ G: ds p a (x,t) £ w iS (x + Ci At,t) Ci , (9) 

where s = 0 or 1 for the fluid and the solid phases, respectively. The interaction 
strength between each fluid and solid surface was adjusted by the parameter G a ads . In 
the simulation, the hydrophilic and hydrophobic regions of the solid particles have 
different adhesion parameter G a ads . The intrinsic contact angle of the hydrophobic 
regions is 1 17.7 degrees, whereas the intrinsic contact angle of the hydrophilic regions 
is 32.8 degrees. 

The macroscopic velocity if ,eq of the o th component in the multiphase model is 
defined by 



-u +- 



(10) 



where u' is the lattice velocity, defined as 



The averaged fluid momentum before and after collision is 2S 



(ii) 



(12) 



where u is the overall fluid velocity, and p = p a is the total mass density of the 
whole fluid. The pressure of the whole fluid is given by 28 

P(x) = c] P a (*) + %p a (x)p°(x). (13) 



The Shan-Chen model modifies the equilibrium velocity in the collision operator to 
include an interactive force, and thus it can realize phase separation and surface- 
tension effects 29 . The wettability of the solid surface (particle and channel) can be 
obtained by adjusting the liquid-solid interaction 30 " 32 . Different contact angles are 
obtained through adjusting the G a ads , as shown in Supplementary Figure S2. The 
evolution computation is parallelized using a Message Passing Interface (MPI). The 
model has been successfully applied in many multiphase investigations, such as 
contact angle 30 33 , contact line dynamics 34 ' 35 , and the effects of wall wettability, topo- 
graphy and micro -structure on drag reduction for fluid flow through micro-chan- 
nels 36 . 
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